!
!  This file is part of REAR (v. 1.0)
!
!  REAR is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, either version 3 of the License, or
!  (at your option) any later version.
!
!  REAR is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with REAR.  If not, see <http://www.gnu.org/licenses/>.
!
! ------------------------------------------------------------------
!
! File "hrm.f90" contains the following program units:
!
!	- Subroutine "PLEgendre":
!              Legendre polynomials (from SHTOOLS) 
!              see: <http://shtools.ipgp.fr>		
!
!       - Subroutine "FINDPX":
!	       Coordinates of the pixels of the Tegmark grid        		
!
!       - subroutine "vector2pixel" (and its dependencies):
!	       Icosahedron Package for Pixelizing the Sphere 
!              written entierely by Mark Tegmark
!              see: <http://space.mit.edu/home/tegmark/icosahedron.html>	
!              
!       - Functions "SINDD" and "COSDD"
!              Real*8 cosine and sine (argument in degres)
!
!       - subroutine "count_rows":
!          Counts the number of rows in a text file
!
!
! Last review by GS on Aug 9, 2013. 
!
! ------------------------------------------------------------------
!
!
!
!
!
!
subroutine PLegendre(p, lmax, z)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!	This subroutine evalutates all of the unnormalized Legendre polynomials 
!	up to degree lmax. 
!
!	Calling Parameters:
!		Out
!			p:	A vector of all unnormalized Legendgre polynomials evaluated at 
!				z up to lmax. The lenght must by greater or equal to (lmax+1).
!		IN
!			lmax:	Maximum degree to compute.
!			z:	Value within [-1, 1], cos(colatitude) or sin(latitude).
!
!	Notes:
!	
!	1.	The integral of Pl**2 over (-1,1) is 2/(2l+1).
!	2.	Values are calculated accoring to the following recursion scheme:
!			P_0(z) = 1.0, P_1(z) = z, and 
!			P_l(z) = (2l-1) * z * P_{l-1}(z) / l - (l-1) * P_{l-2}(z) / l
!
!	Dependencies:	None
!
!	Written by Mark Wieczorek June 2004
!
! ----> Modified to SINGLE PRECISION by Giorgio Spada 2007 
! ----> Also modified for the management of ERROR conditions 
!
!	Copyright (c) 2005, Mark A. Wieczorek
!	All rights reserved.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	
	implicit none
	integer, intent(in) ::	lmax
	real*8 paux(lmax+1)	
	real*8, intent(out) ::	p(0:lmax+1)		
       	real*8, intent(in) ::	z
       	real*8 :: pm2, pm1, pl
      	integer :: l, j

	if(size(p) < lmax+1) then
 	        Write(* ,*) "Error --- PlegendreL" 
		Write(* ,*) "P must be dimensioned as (LMAX+1) where LMAX is ", lmax
		Write(* ,*) "Input array is dimensioned ", size(p)
		Write(* ,*) "The program will STOP ----------------"	   
	        Stop
     	elseif (lmax < 0) then 
 	        Write(* ,*) "Error --- PlegendreL" 
		Write(* ,*) "LMAX must be greater than or equal to 0."
		Write(* ,*) "Input value is ", lmax
		Write(* ,*) "The program will STOP ----------------"	   
     	elseif(abs(z) > 1.) then
 	        Write(* ,*) "Error --- PlegendreL" 
		Write(* ,*) "ABS(Z) must be less than or equal to 1."
		Write(* ,*) "Input value is ", z
		Write(* ,*) "The program will STOP ----------------"	   
	        Stop
     	endif
      	
   	pm2  = 1.
      	paux(1) = 1.
      	
      	pm1  = z
      	paux(2) = pm1
      	
      	do l = 2, lmax
         	pl = (float(2*l-1) * z * pm1 - float(l-1) * pm2)/float(l)
         	paux(l+1) = pl
         	pm2  = pm1
         	pm1  = pl
      	enddo
!
        do j=0, lmax
             p(j)=paux(j+1)
	enddo
!
end subroutine PLegendre
!
!
!
!
!
!
	Subroutine FINDPX(RES,NP,LONP,LATP)
!
	implicit NONE 
!
! # Given resolution RES, this routine returns the number of pixels NP,
!   and the longitude and latitude arrays LONP, LATP (1:NP) with the pixels
!   coordinates. THis done through a call to "pixel2vector" by M. Tegmark- 
!   
!   Written by GS on October 2, 2008 for the implementation of the 
!   3D velocity maps for SELEN 2.7 - but perhaps useful also for other 
!   purposes... 
!
	INTEGER J, NP, RES
        REAL*4 R(0:19,3,3), V(0:11,3), VECT(3), X,Y,Z  
	REAL*4 LONP(1:NP), LATP(1:NP), LONPX, LATPX  
	REAL*4,  PARAMETER :: PI=3.14159265358979323840 
!
        call compute_matrices (r)
        call compute_corners  (v)
!
	np=40*res*(res-1)+12 
!
        do 1 j=0,np-1
	  	call pixel2vector (j,res,r,v,vect)
	        x=vect(1)
		y=vect(2)
		z=vect(3)
! 
! --- Polar pixel 
          	if(x==0..and.y==0.) then 
          		lonpx=0.
                	if(z>=0.) latpx = +90.
                	if(z<=0.) latpx = -90.
! --- Ordinary pixel 
	  		else                         	  
	  		lonpx = atan2(y,x)*180./pi 
          		if (lonpx < 0.) lonpx=360.+ lonpx
          		latpx = 90.-acos(z)*180./pi 
          	endif
	 lonp(j+1)=lonpx 
	 latp(j+1)=latpx 	
!
1       continue 
!
	end subroutine FINDPX 
!
!
!
!
!
!
!
!								      
!				   ~~~ 				      
!				 --- --- 		              
!	 		-+-+-+-+-+-+-+-+-+-+-+-+-		      
! ----===========================================================---- 
!         The code which follows is entierely by Mark Tegmark ...     
! ----===========================================================---- 
! 			-+-+-+-+-+-+-+-+-+-+-+-+-		      
!				 --- --- 			      
!				   ~~~				      



!	
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!     THESE SUBROUTINES ARE ALL YOU NEED TO CALL FROM OUTSIDE    !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	
!!!	These subroutines convert between unit vectors and         !!!
!!!     pixel numbers, and are the only ones that the user of this !!!
!!!     package calls repeatedly:				   !!!
!!!	  subroutine vector2pixel(vector,resolution,R,v,pixel)     !!!
!!!	  subroutine pixel2vector(pixel,resolution,R,v,vector)     !!!
!!!                                                                !!!
!!!	These subroutines are called only once, in the beginning,  !!!
!!!     and compute the necessary rotation matrices and corner     !!!
!!!     vectors once and for all:                                  !!!
!!!	  subroutine compute_matrices(R)                           !!!
!!!	  subroutine compute_corners(v)                            !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	

	subroutine vector2pixel(vector,resolution,R,v,pixel)
	implicit none
	real vector(3), R(0:19,3,3), v(0:11,3)
	real A(3,3), vec(3), x, y
	integer resolution, pixel
	integer pix, face, pixperface, ifail
	if (resolution.lt.1) then
		write(*,*) "ERROR in the Tegmark modules: Resolution must exceed 0"
		STOP
	Endif
	pixperface = 2*resolution*(resolution-1)
	call find_face(vector,R,face)
	call getmatrix(face,R,A)
	call vecmatmul2(A,vector,vec)	
	x	= vec(1)/vec(3)
	y	= vec(2)/vec(3)
	call adjust(x,y)
	call tangentplanepixel(resolution,x,y,pix,ifail)
	if (ifail.gt.0) then 
	  ! Try the runner-up face:
	  call find_another_face(vector,R,face)
	  call getmatrix(face,R,A)
	  call vecmatmul2(A,vector,vec)	
	  x	= vec(1)/vec(3)
	  y	= vec(2)/vec(3)
	  call adjust(x,y)
	  call tangentplanepixel(resolution,x,y,pix,ifail)
	end if	
	pixel = face*pixperface + pix
	if (ifail.gt.0) then
	  ! The pixel wasn't in any of those two faces,
	  ! so it must be a corner pixel.
	  call find_corner(vector,v,pix)
	  pixel = 20*pixperface + pix
	end if
	return
	end

	subroutine pixel2vector(pixel,resolution,R,v,vector)
	! Returns a unit vector pointing towards pixel.
	! Resolution must be an even, positive integer.
	implicit none
	real vector(3), R(0:19,3,3), v(0:11,3)
	real A(3,3), x, y, norm
	integer resolution, pixel
	integer pix, face, pixperface
	if (resolution.lt.1) then 
		Write(*,*) "ERROR in the Tegmark modules: Resolution must exceed 0"
		STOP
	ENDIF
	pixperface = 2*resolution*(resolution-1)
	if (pixel.lt.0) then 
		Write(*,*)"ERROR in the Tegmark modules: negative pixel number"	
		STOP
	ENDIF
	if (pixel.ge.20*pixperface+12) then 
		Write(*,*)"ERROR in the Tegmark modules: pixel number too large"
		STOP 
	Endif
	if (pixperface.gt.0) then 
	  face = pixel/pixperface
	  if (face.gt.20) face = 20
	else ! There are no pixels at all on the faces - it's all just corners.
	  face = 20
	end if	
	pix = pixel - face*pixperface
	if (face.lt.20) then
	  ! The pixel is on one of the 20 faces:
	  call tangentplanevector(pix,resolution,x,y)
	  call unadjust(x,y)
	  norm 	= sqrt(x*x+y*y+1)
	  vector(1)	= x/norm
	  vector(2)	= y/norm
	  vector(3)	= 1./norm
	  call getmatrix(face,R,A)
	  call vecmatmul1(A,vector,vector)
	else
	  ! This is a corner pixel:
	  if (pix.gt.11) then 
	  	Write(*,*)"ERROR in the Tegmark modules: pixel number too big"
	  	STOP 
	  endif
	  vector(1) = v(pix,1)
	  vector(2) = v(pix,2)
	  vector(3) = v(pix,3)
	end if
	return
	end

	subroutine compute_matrices(R)
	! On exit, R will contain the 20 rotation matrices
	! that rotate the 20 icosahedron faces 
	! into the tangent plane
	! (the horizontal plane with z=1). 
	! Only called once, so speed is irrelevant.
	implicit none
	real R(0:19,3,3)
	real A(3,3), B(3,3), C(3,3), D(3,3), E(3,3)
	real pi, sn, cs, ct, x
	integer i,j,n
	do i=1,3
	  do j=1,3
	    A(i,j) = 0.
	    B(i,j) = 0.
	    C(i,j) = 0.
	    D(i,j) = 0.
	  end do
	end do
	pi 	= 4.*atan(1.)
	x 	= 2.*pi/5.
	cs 	= cos(x)
	sn	= sin(x)
	A(1,1)	= cs
	A(1,2)	= -sn
	A(2,1)	= sn
	A(2,2)	= cs
	A(3,3) 	= 1.
	! A rotates by 72 degrees around the z-axis.
	x 		= pi/5.
	ct 		= cos(x)/sin(x)
	cs		= ct/sqrt(3.)
	sn		= sqrt(1-ct*ct/3.)
	C(1,1)	= 1
	C(2,2)	= cs
	C(2,3)	= -sn
	C(3,2)	= sn
	C(3,3) 	= cs
	! C rotates around the x-axis so that the north pole 	
	! ends up at the center of face 1.
	cs		= -0.5
	sn		= sqrt(3.)/2
	D(1,1)	= cs
	D(1,2)	= -sn
	D(2,1)	= sn
	D(2,2)	= cs
	D(3,3)	= 1.
	! D rotates by 120 degrees around z-axis.
	call matmul1(C,D,E)
	call matmul2(E,C,B)	! B = CDC^t
	! B rotates face 1 by 120 degrees. 
	do i=1,3
	  do j=1,3
	    E(i,j) = 0.
	  end do
	  E(i,i) = 1.
	end do	! Now E is the identity matrix.
	call putmatrix(0,R,E)	
	call matmul1(B,A,E)
	call matmul1(B,E,E)
	call putmatrix(5,R,E)	
	call matmul1(E,A,E)
	call putmatrix(10,R,E)
	call matmul1(E,B,E)
	call matmul1(E,B,E)
	call matmul1(E,A,E)
	call putmatrix(15,R,E)
	do n=0,15,5	
	  call getmatrix(n,R,E)
	  do i=1,4
	    call matmul1(A,E,E)
	    call putmatrix(n+i,R,E)
	  end do
	end do
	! Now the nth matrix in R will rotate 
	! face 1 into face n. 
	! Multiply by C so that they will rotate
	! the tangent plane into face n instead:
	do n=0,19
	  call getmatrix(n,R,E)
	  call matmul1(E,C,E)
	  call putmatrix(n,R,E)
	end do
	return
	end

	subroutine compute_corners(v)
	! On exit, v will contain unit vectors pointing toward 
	! the 12 icoshedron corners. 
	implicit none
	real v(0:11,3)
	real pi, z, rho, dphi
	integer i
	pi = 4.*atan(1.)
	dphi = 2.*pi/5.
	! First corner is at north pole:
	v(0,1) = 0.
	v(0,2) = 0.
	v(0,3) = 1.
	! The next five lie on a circle, with one on the y-axis:
	z = 0.447213595	! This is 1/(2 sin^2(pi/5)) - 1
	rho = sqrt(1.-z*z)
	do i=0,4
	  v(1+i,1) = -rho*sin(i*dphi)
	  v(1+i,2) =  rho*cos(i*dphi)
	  v(1+i,3) = z 
	end do
	! The 2nd half are simply opposite the first half:
	do i=0,5
	  v(6+i,1) = -v(i,1)
	  v(6+i,2) = -v(i,2)
	  v(6+i,3) = -v(i,3)
	end do
	return
	end

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!     THE SUBROUTINES BELOW ARE SUBORDINATE TO THOSE ABOVE, AND  !!!
!!!     CAN BE SAFELY IGNORED BY THE GENERAL USER OF THE PACKAGE.  !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	
!!!	These subroutines perform some standard linear algebra:    !!!
!!!	  subroutine matmul1(A,B,C)                                !!!
!!!	  subroutine matmul2(A,B,C)                                !!!	
!!!	  subroutine matmul3(A,B,C)                                !!!	
!!!	  subroutine vecmatmul1(A,b,c)	                           !!!	
!!!	  subroutine vecmatmul2(A,b,c)	                           !!!	
!!!                                                                !!!
!!!	These subroutines copy matrices in and out of storage:     !!!
!!!	  subroutine getmatrix(n,R,A)                              !!!
!!!	  subroutine putmatrix(n,R,A)                              !!!
!!!                                                                !!!
!!!     These subroutines help vector2pixel reduce the 3D sphere   !!!
!!!     problem to a problem on an equilateral triangle in the     !!!
!!!     z=1 tangent plane (an icosahedron face):                   !!!		
!!!	  subroutine find_face(vector,R,face)                      !!!
!!!	  subroutine find_another_face(vector,R,face)              !!!
!!!	  subroutine find_corner(vector,v,corner)                  !!!
!!!                                                                !!!
!!!     These subroutines pixelize this triangle with a regular    !!!
!!!     triangular grid:                                           !!!
!!!	  subroutine find_mn(pixel,resolution,m,n)                 !!!
!!!	  subroutine tangentplanepixel(resolution,x,y,pix,ifail)   !!!
!!!	  subroutine tangentplanevector(pix,resolution,x,y)        !!!
!!!                                                                !!!
!!!     These subroutines reduce the area equalization problem to  !!!
!!!     one on the right triangle in the lower right corner:       !!!
!!!	  subroutine find_sixth(x,y,rot,flip)                      !!!
!!!	  subroutine rotate_and_flip(rot,flip,x,y)	           !!!
!!!	  subroutine adjust(x,y)                                   !!!
!!!	  subroutine unadjust(x,y)                                 !!!
!!!                                                                !!!
!!!     These subroutines perform the area equalization mappings   !!!
!!!     on this right triangle:                                    !!!
!!!	  subroutine adjust_sixth(x,y)                             !!!
!!!	  subroutine unadjust_sixth(x,y)                           !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	

	subroutine matmul1(A,B,C)	
	! Matrix multiplication C = AB.
	! A, B and C are allowed to be physiclly the same.
	implicit none
	real A(3,3), B(3,3), C(3,3), D(3,3), sum
	integer i,j,k
	sum = 0.
	do i=1,3
	  do j=1,3
	    sum = 0.
	    do k=1,3
	      sum = sum + A(i,k)*B(k,j)
	    end do
	    D(i,j) = sum
	  end do
	end do
	call copymatrix(D,C)
	return
	end

	subroutine matmul2(A,B,C)	
	! Matrix multiplication C = AB^t
	! A, B and C are allowed to be physically the same.
	implicit none
	real A(3,3), B(3,3), C(3,3), D(3,3), sum
	integer i,j,k
	sum = 0.
	do i=1,3
	  do j=1,3
	    sum = 0.
	    do k=1,3
	      sum = sum + A(i,k)*B(j,k)
	    end do
	    D(i,j) = sum
	  end do
	end do
	call copymatrix(D,C)
	return
	end

	subroutine matmul3(A,B,C)	
	! Matrix multiplication C = A^t B
	! A, B and C are allowed to be physically the same.
	implicit none
	real A(3,3), B(3,3), C(3,3), D(3,3), sum
	integer i,j,k
	sum = 0.
	do i=1,3
	  do j=1,3
	    sum = 0.
	    do k=1,3
	      sum = sum + A(k,i)*B(k,j)
	    end do
	    D(i,j) = sum
	  end do
	end do
	call copymatrix(D,C)
	return
	end

	subroutine vecmatmul1(A,b,c)	
	! Matrix multiplication c = Ab
	! b and c are allowed to be physically the same.
	implicit none
	real A(3,3), b(3), c(3), d(3), sum
	integer i,j
	sum = 0.
	do i=1,3
	  sum = 0.
	  do j=1,3
	    sum = sum + A(i,j)*b(j)
	  end do
	  d(i) = sum
	end do
	call copyvector(d,c)
	return
	end

	subroutine vecmatmul2(A,b,c)	
	! Matrix multiplication c = A^tb
	! b and c are allowed to be physiclly the same.
	implicit none
	real A(3,3), b(3), c(3), d(3), sum
	integer i,j
	sum = 0.
	do i=1,3
	  sum = 0.
	  do j=1,3
	    sum = sum + A(j,i)*b(j)
	  end do
	  d(i) = sum
	end do
	call copyvector(d,c)
	return
	end

	subroutine copymatrix(A,B)	
	! B = A
	implicit none
	real A(3,3), B(3,3)
	integer i,j
	do i=1,3
	  do j=1,3
	    B(i,j) = A(i,j)
	  end do
	end do
	return
	end

	subroutine copyvector(a,b)	
	! b = a
	implicit none
	real a(3), b(3)
	integer i
	do i=1,3
	  b(i) = a(i)
	end do
	return
	end

	subroutine getmatrix(n,R,A)	
	! A = the nth matrix in R
	implicit none
	real R(0:19,3,3), A(3,3)
	integer i,j,n
	do i=1,3
	  do j=1,3
	    A(i,j) = R(n,i,j)
	  end do
	end do
	return
	end

	subroutine putmatrix(n,R,A)	
	! the nth matrix in R = A
	implicit none
	real R(0:19,3,3), A(3,3)
	integer i,j,n
	do i=1,3
	  do j=1,3
	    R(n,i,j) = A(i,j)
	  end do
	end do
	return
	end

	subroutine find_face(vector,R,face)
	! Locates the face to which vector points.
	! Computes the dot product with the vectors
	! pointing to the center of each face and picks the
	! largest one. 
	! This simple routine can be substantially accelerated
	! by adding a bunch of if-statements, to avoid looping
	! over more than a few faces.
	implicit none
	real vector(3), R(0:19,3,3), dot, max
	integer n,face,i
	max = -17.
	do n=0,19
	  dot = 0.
	  do i=1,3
	    dot = dot + R(n,i,3)*vector(i)
	  end do
	  if (dot.gt.max) then
	    face = n
	    max = dot
	  end if
	end do
	return
	end

	subroutine find_another_face(vector,R,face)
	! Computes the dot product with the vectors
	! pointing to the center of each face and picks the
	! largest one other than face.
	! This simple routine can be substantially accelerated
	! by adding a bunch of if-statements, to avoid looping
	! over more than a few faces.
	implicit none
	real vector(3), R(0:19,3,3), dot, max
	integer n,face,facetoavoid,i
	facetoavoid = face
	max = -17.
	do n=0,19
	  if (n.ne.facetoavoid) then 
	    dot = 0.
	    do i=1,3
	      dot = dot + R(n,i,3)*vector(i)
	    end do
	    if (dot.gt.max) then
	      face = n
	      max = dot
	    end if
	  end if
	end do
	return
	end

	subroutine find_corner(vector,v,corner)
	! Locates the corner to which vector points.
	! Computes the dot product with the vectors
	! pointing to each corner and picks the
	! largest one. 
	! This simple routine can be substantially accelerated
	! by adding a bunch of if-statements, but that's pretty
	! pointless since it gets called so rarely.
	implicit none
	real vector(3), v(0:11,3), dot, max
	integer corner,n,i
	max = -17.
	do n=0,11
	  dot = 0.
	  do i=1,3
	    dot = dot + v(n,i)*vector(i)
	  end do
	  if (dot.gt.max) then
	    corner = n
	    max = dot
	  end if
	end do
	return
	end

	subroutine find_mn(pixel,resolution,m,n)
	! Computes the integer coordinates (m,n) of the pixel 
	! numbered pix on the basic triangle.
	implicit none
	integer pixel, pix, resolution, m, n
	integer interiorpix , pixperedge
	pix 	    = pixel
	interiorpix = (2*resolution-3)*(resolution-1)
	pixperedge  = (resolution)-1
	if (pix.lt.interiorpix) then 
	  ! The pixel lies in the interior of the triangle.
	  m = (sqrt(1.+8.*pix)-1.)/2. + 0.5/resolution
	  ! 0.5/resolution was added to avoid problems with
	  ! rounding errors for the case when n=0.
	  ! As long as you don't add more than 2/m, you're OK.
	  n = pix - m*(m+1)/2
	  m = m + 2
	  n = n + 1
	  goto 555
	end if
	pix = pix - interiorpix 
	if (pix.lt.pixperedge) then
	  ! The pixel lies on the bottom edge.
	  m = 2*resolution-1
	  n = pix+1
	  goto 555
	end if
	pix = pix - pixperedge
	if (pix.lt.pixperedge) then
	  ! The pixel lies on the right edge.
	  m = 2*resolution-(pix+2)
	  n = m
	  goto 555
	end if
	pix = pix - pixperedge
	! The pixel lies on the left edge.
	m = pix+1
	n = 0
555	return
	end

	subroutine tangentplanepixel(resolution,x,y,pix,ifail)
	! Finds the hexagon in which the point (x,y) lies
	! and computes the corresponding pixel number pix.
	! Returns ifail=0 if (x,y) lies on the face, 
	! otherwise returns ifail=1.
	implicit none
	real x, y, a, b, c, d, edgelength
	parameter (c=0.866025404)	! sqrt(3)/2
	parameter(edgelength=1.3231690765)
	! The edge length of the icosahedron is 
	! sqrt(9 tan^2(pi/5) - 3) when scaled so that
	! it circumscribes the unit sphere. 
	integer resolution, pix, ifail, i, j, k, m, n, r2
	r2	= 2*resolution
	a 	= 0.5*x
	b 	= c*y
	d 	= 0.5*edgelength/r2
	i 	= x/d 	  + r2
	j 	= (a+b)/d + r2
	k 	= (a-b)/d + r2
	m 	= (r2+r2-j+k-1)/3
	n 	= (i+k+1-r2)/3
	pix 	= (m-2)*(m-1)/2 + (n-1)
	ifail = 0
	if (m.eq.r2-1) then		! On bottom row
	  if ((n.le.0).or.(n.ge.resolution)) then
	    ifail=1
	  end if
	  goto 666				! Pix already correct
	end if
	if (n.eq.m) then			! On right edge
	  k = (r2-1) - m
	  if ((k.le.0).or.(k.ge.resolution)) then
	    ifail = 1
	  else
	    pix = (r2-2)*(resolution-1) + k - 1
 	  end if
	  goto 666
	end if
	if (n.eq.0) then			! On left edge
	  if ((m.le.0).or.(m.ge.resolution)) then
	    ifail = 1
	  else
	    pix = (r2-1)*(resolution-1) + m - 1
	  end if
	end if
666	return
	end

	subroutine tangentplanevector(pix,resolution,x,y)
	! Computes the coordinates (x,y) of the pixel 
	! numbered pix on the basic triangle.
	implicit none
	real x, y, c1, c2, edgelength
	parameter(c1=0.577350269)	! 1/sqrt(3)
	parameter(c2=0.866025404)	! sqrt(3)/2
	parameter(edgelength=1.3231690765)
	! The edge length of the icosahedron is 
	! sqrt(9 tan^2(pi/5) - 3) when scaled so that
	! it circumscribes the unit sphere. 
	integer pix, resolution, m, n
	call find_mn(pix,resolution,m,n)
	x	= edgelength*(n-0.5*m)/(2*resolution-1)
	y 	= edgelength*(c1-(c2/(2*resolution-1))*m)
	return
	end

	subroutine find_sixth(x,y,rot,flip)
	! Find out in which sixth of the basic triangle
	! the point (x,y) lies, identified by the
	! two integers rot (=0, 1 or 2) and flip = 0 or 1).
	! rot and flip are defined such that the sixth is 
	! mapped onto the one at the bottom right by
	! these two steps:
	! 1. Rotate by 120 degrees anti-clockwise, rot times.
	! 2. Flip the sign of x if flip = 1, not if flip=0.
	! The if-statements below go through the six cases 
	! anti-clockwise, starting at the bottom right.
	implicit none
	real x, y, c, d
	parameter(c=1.73205081)		! sqrt(3)
	integer rot, flip
	d = c*y
	if (x.ge.0) then
	  if (x.le.-d) then
	    rot  = 0
	    flip = 0
	  else
	    if (x.ge.d) then
	      rot  = 2
	      flip = 1
	    else 
	      rot  = 2
	      flip = 0
	    end if
	  end if
	else
	  if (x.ge.-d) then
	    rot  = 1
	    flip = 1
	  else
	    if (x.le.d) then
	      rot  = 1
	      flip = 0
	    else 
	      rot  = 0
	      flip = 1
	    end if
	  end if
	end if
	return
	end

	subroutine rotate_and_flip(rot,flip,x,y)
	implicit none
	real x, y, x1, cs, sn, c
	integer rot, flip
	parameter(cs=-0.5)
	parameter(c=0.866025404)	! sqrt(3)/2
	if (rot.gt.0) then
	  if (rot.eq.1) then
	    sn = c	! Rotate 120 degrees anti-clockwise 
	  else 
	    sn = -c	! Rotate 120 degrees anti-clockwise 
	  end if
	  x1 = x
	  x  = cs*x1 - sn*y
	  y  = sn*x1 + cs*y
	end if
	if (flip.gt.0) x = -x
	return
	end

	subroutine adjust(x,y)
	! Maps the basic triangle onto itself in such a way
	! that pixels will have equal area when mapped onto
	! the sphere. 
	implicit none
	real x, y
	integer rot, flip
	call find_sixth(x,y,rot,flip)
	call rotate_and_flip(rot,flip,x,y)
	call adjust_sixth(x,y)
	! Now rotate & flip the sixth back into its
	! original position:
	if ((flip.eq.0).and.(rot.gt.0)) then  
	  call rotate_and_flip(3-rot,flip,x,y)
	else
	  call rotate_and_flip(rot,flip,x,y)
	end if
	return
	end
	
	subroutine unadjust(x,y)
	! Performs the inverse of what adjust does. 
	implicit none
	real x, y
	integer rot, flip
	call find_sixth(x,y,rot,flip)
	call rotate_and_flip(rot,flip,x,y)
	call unadjust_sixth(x,y)
	! Now rotate & flip the sixth back into its
	! original position:
	if ((flip.eq.0).and.(rot.gt.0)) then  
	  call rotate_and_flip(3-rot,flip,x,y)
	else
	  call rotate_and_flip(rot,flip,x,y)
	end if
	return
	end
	
	subroutine adjust_sixth(x,y)
	! Maps the basic right triangle (the sixth of the face that
	! is in the lower right corner) onto itself in such a way
	! that pixels will have equal area when mapped onto the sphere. 
	implicit none
	real x, y, u, v, g, v2, root, trig, eps, scale
	parameter(eps=1.e-14, scale=1.09844)
	parameter(g=1.7320508075689)		! sqrt(3)
	u 		= x  + eps
	v		= -y + eps
	v2		= v*v
	root		= sqrt(1.+4.*v2)
	trig		= atan((g*root-g)/(root+3.))
	y		= sqrt(trig*2./g)
	x		= sqrt((1.+4.*v2)/(1.+u*u+v2))*u*y/v
	x		= scale*x
	y 		= -scale*y	
	return
	end
	
	subroutine unadjust_sixth(x,y)
	! Performs the inverse of what adjust_sixth does.
	implicit none
	real x, y, u, v, g, v2, y2, tmp, trig, eps, scale
	parameter(eps=1.e-14, scale=1.09844)
	parameter(g=1.7320508075689)		! sqrt(3)
	u 		=  x/scale + eps
	v		= -y/scale + eps
	v2		= v*v
	trig		= tan(g*v2/2.)
	tmp		= (g+3.*trig)/(g-trig)
	y2		= (tmp*tmp-1.)/4.
	y		= sqrt(y2)
	tmp		= v2*(1.+4.*y2) - u*u*y2
	x	 	= u*y*sqrt((1.+y2)/tmp) 
	y 		= -y	
	return
	end

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!     END OF THE ICOSAHEDRON PACKAGE FOR PIXELIZING THE SPHERE   !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	
!
!
!
!
!
!
!
!
	FUNCTION SINDD(ALFA)
	implicit NONE
!
! --- Assuming that ALFA is given in DEGREES, it returns the SINE of ALFA
!     *********************** GS and FC 11-02-2009 **********************
!	
	REAL*8,  PARAMETER :: PI=3.14159265358979323840 
	REAL*8 SINDD, ALFA, ALFARAD              
	ALFARAD=ALFA*PI/180.	
	SINDD=SIN(ALFARAD)
	
	END FUNCTION SINDD 
!
!
!
!
!
!
	FUNCTION COSDD(ALFA)
	implicit NONE
!
! --- Assuming that ALFA is given in DEGREES, it returns the COSINE of ALFA
!     ************************ GS and FC 11-02-2009 ***********************
!	
	REAL*8,  PARAMETER :: PI=3.14159265358979323840 
	REAL*8 COSDD, ALFA, ALFARAD               
	ALFARAD=ALFA*PI/180.	
	COSDD=COS(ALFARAD)	
!
	END FUNCTION COSDD
!
!
!
!
!
    subroutine count_rows(lun,n)
    implicit none
!
! --- Returns as N the number of rows in text file opened at LUN
!     DM 28.04.2014
!
    integer lun,n
!
    rewind(lun)
!    
    n=0
901 read(lun,*,end=902)
    n=n+1
    go to 901
    902 continue
!
    return
!
    end subroutine count_rows
!    
!
!
